How to Calculate Tortuosity Easily? 
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Abstract. Tortuosity is one of tiie key parameters describing the geometry and transport properties of porous media. It is 
defined eitlier as an average elongation of fluid patfis or as a retardation factor tfiat measures tlie resistance of a porous 
medium to tlie flow. However, in contrast to a retardation factor, an average fluid patfi elongation is difficult to compute 
numerically and, in general, is not measurable directly in experiments. We review some recent achievements in bridging the 
gap between the two formulations of tortuosity and discuss possible method of numerical and an experimental measurements 
of the tortuosity directly from the fluid velocity field. 
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1. INTRODUCTION 

One of the main problems in porous media physics is to 
find out how the value of permeability, which synthet- 
ically describes flow retardation by the porous medium 
structure, is correlated with the geometry of the medium. 
Another problem is to define macroscopic parameters 
that could be used to distinguish various kinds of porous 
media. One of the parameters that is often used for both 
of these purposes is the tortuosity. 

The notion of tortuosity was introduced to porous 
media studies by Carman [1], who considered a flow 
through a bed of sand and proposed the tortuosity as 
a factor that accounts for effective elongation of fluid 
paths. Assuming that a porous bed of thickness L can be 
regarded as a bundle of capillaries of equal length Leff 
and constant cross-section, he proposed the following 
semi-empirical Kozeny-Carman formula [1,2, 3] 
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which relates the permeability {k) to four structural pa- 
rameters: the porosity (p, the specific surface area S, the 
shape factor /3, and the hydraulic tortuosity T, 



FIGURE 1. Streamlines in the fluid flow through three- 
dimensional random model of porous media at porosity (p = 0.6 
and tortuosity T=1.15. 



porous medium filled with this fluid. Within the simple 
capillary model, T^\ is related to L^ft/L through 
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The simple capillary model used by Kozeny and Car- 
man can be easily applied to other forms of transport 
through porous media. For example, the electric tortu- 
osity (Tel) is defined as a retardation factor [4], 
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where Cfi is the electrical conductivity of a conductive 
fluid and (7p is the effective electrical conductivities of a 



and a similar relation holds for the diffusive tortuosity 
[4]. 

Comparison of Eqs. (2) and (4), as well as a research 
into the literature, reveal the first problem with tortuosity: 
depending on the context, this term can be related to T, 
T^ or even T^' or T^^ [2, 3, 4, 5]. The second problem 
is that a link between the tortuosity defined as an average 
elongation of fluid paths, as in Eq. (2) (see Fig. 1), and 
the tortuosity defined as a retardation factor, as in Eq. (4), 
is well-defined only for the capillary model. It is not clear 
that a similar correlation exists for arbitrary porous me- 



dia. The third problem is an imprecise definition of the 
effective fluid path length (Leg) in Eq. (2). In real porous 
media flow paths are extremely complicated, as the fluid 
fluxes continuously change in sectional area, shape and 
orientation as well as branch and rejoin, and this observa- 
tion led several researchers [3] to believe that Leff can be 
defined only in relatively simple network models, which, 
however, can be analyzed without this notion. The fourth 
problem is that the Carman-Kozeny formula, Eq. (1), ac- 
tually defines the product j3r^ (known as the Kozeny 
constant) and if the tortuosity could not be measured for 
general porous media, the shape factor and hydraulic tor- 
tuosity would become essentially indeterminate quanti- 
ties, rendering the tortuosity a 'fudge factor' used to fit 
the model to the experimental data [3, 6]. 

In this paper we discuss some recent achievements in 
solving the above-mentioned problems. In particular, we 
show some applications and impUcations of our recent 
formula for the tortuosity [7] 
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where (u) is the average magnitude of the intrinsic ve- 
locity over the entire system volume and (ux) is the vol- 
umetric average of its component along the macroscopic 
flow direction. 



2. SOLUTION TO THE FLUID FLOW 
PROBLEM 

In order to compute the tortuosity defined with Eq. (5) 
one has to know the steady state velocity field. This 
may be accomplished either experimentally, e.g. by using 
the particle image velocimetry methods [8, 9, 10], the 
magnetic resonanse imaging [11, 12] or numerically, by 
finding the solution to the Navier-Stokes equations in the 
pore space of a porous medium. 

Here we take the numerical approach. We use the 
lattice Boltzmann method (LBM) [13]. In this method, 
which originates from the kinetic theory of gases and the 
lattice gas automata models [14], the fluid is modeled 
as consisting of fictive particles propagating and collid- 
ing on a discrete lattice. Space, time and velocities are 
all discrete, with space usually discretized into a regular 
grid, time discretized into equal intervals, and velocities 
restricted to just a few vectors C; related to the geometry 
of the lattice. The state of the system is fully character- 
ized by distribution functions /i(x,r) g [0, 1] defined for 
each lattice node x, discrete time f, and c,. They can be 
interpreted as being proportional to the number of par- 
ticles that at time f are at node x and have velocity c,. 
It was shown that solving the LBM model leads to the 



solution of the incompressible Navier-Stokes equations 

[15]. 



2.1. Sailfish and tortuosity computation 

All our numerical computations were performed using 
the Sailfish library [16], which is an implementation of 
the LBM method running on graphics processing units 
(GPUs), an emerging platform for high-performance 
computing. Sailfish is an open-source project written in 
the Python programming language, and hence is a highly 
customisable software. In particular, implementation of 
Eq. (5) in Sailfish is trivial and consists of just a few lines 
of code. 



3. TORTUOSITY COMPUTATION 

A general discussion of Eq. (5), its relation to Eq. (2) 
and conditions of applicability can be found in [7]. Here 
we use the fact that the LBM method uses a regular (i.e. 
cubic) mesh. This allows to approximate Eq. (5) with 
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where r runs through all lattice nodes. Note that this 
simple formula can be used not only in numerical studies, 
but also for the data obtained experimentally. 

In subsequent subsections we test how Eq. (6) behaves 
for flows in various geometries of increasing complexity. 



3.1. Inclined channel 

The simplest model of a porous medium approximates 
it as a bunch of straight pipes, each inclined at an angle 
a to the macroscopic flow direction. In this case all 
streamlines are of the same length and Eq. (2) yields 
T = 1/cosa. We constructed two two-dimensional (2D) 
configurations of an inclined channel of height h, with 
/i = 5 and 20 lattice units (l.u.) (see Fig. 2). We used the 
mesh resolution 240 x 800 (l.u.), set the lattice kinematic 
viscosity at v = 0.01 and assumed the lattice velocity 
Mmax = 0.02 as the maximum value of the developed 
velocity profile at the inlet and outlet. Starting from 
a = 0, we succesively rotated the channel preserving its 
width. For each inclination angle the steady state was 
found using the LBM method and then the tortuosity 
was calculated from Eq. (6). The results are depicted in 
Fig. 3. They agree well with the values expected from 
geometric considerations even for a narrow channel of 
height h=5 l.u. (the relative difference does not exceed 
5%). Small errors observed in our simulations stem from 
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FIGURE 2. A channel rotated by an angle a. 



FIGURE 4. The U-shaped channel with a step depth b. 
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FIGURE 3. Hydraulic tortuosity T, Eq. (5), in a channel of 
height h (lattice units) rotated by an angle a. The solid line 
shows the geometric tortuosity 1/cosa. 




FIGURE 5. Hydraulic tortuosity, T, computed using Eq. (5) 
in a U-shaped channel as a function of the relative step depth 
b/H for two channel heights h (lattice units). The solid line 
represents the geometric tortuosity for the same system. 



discretization errors and the staircased approximation of 
the channel boundaries in the LBM method [17], which 
results in T being slightly overestimated. These eiTors 
decrease with an increasing value of h. 



3.2. U-shaped channel 

Next, we constructed a U-shaped channel geometry 
(see Fig. 4). The mesh resolution was set at L x // = 
300 X 200 (I.U.), and we assumed v = 0.01 and Wmax = 
0.01. We started from the step depth b = Q and increased 
it until b = H /2. For each b two channel heights, /j = 5 
and h ~ 20 (l.u.) were investigated. In this geometry the 
geometric tortuosity Tg = {L + 2b)/L. Comparison of the 
tortuosity determined from Eq. (6) with Tg for various 
values ofb and h is shown in Fig. 5. As expected, a linear 
dependency of geometric and hydraulic tortuosities on b 
is visible, but a deviation of T from T^ at larger chan- 
nel depths b is also noticeable. To understand this effect 
we analyzed the flow streamlines for this systems (data 



not shown) and observed that they follow the geometry 
nicely only at the straight parts of the channel. At cor- 
ners, however, a tendency of flow paths to seek a shorter 
path is visible. This leads to the hydraulic toituosity T 
being smaller than the geometric one, as shown in Fig. 5. 
This effect is more visible for larger h, as in this case the 
corner deformation is more pronounced. 



3.3. Two-dimensional overlapping circles 

The next model we considered was a 2D system in 
which the porous matrix was modeled by circles that 
were free to overlap. We used a 900 x 600 mesh at which 
we randomly deposited monosized circles of radius r = 
10 l.u. No-slip boundary conditions were imposed on the 
top and bottom walls and periodic boundary conditions 
were assumed at the inlet and outlet. The flow was driven 
by an external force field whose magnitude was chosen 
so that the Reynolds number Re < 1 (creeping flow). 




FIGURE 6. Flow streamlines in a realization of a two- 
dimensional model of a porous medium built of overlapping 
circles of radius R = 10 (l.u.) at porosity (p = 0.85. The tortu- 
osity computed using Eq. (6) is T = 1.13. 



The circles were deposited only in the central, 600 x 600 
(l.u.) area of the mesh, see Fig. 6. The remaining space 
was kept empty to minimize the influence of inlet and 
outlet boundary conditions. We ran the simulation for 
40 000 steps until the steady state was reached and used 
Eq. (6) to compute the tortuosity. As seen in Fig. 6, the 
flow streamlines in this model can be quite complex and 
"tortuous". 

In our previous studies [18, 19] we considered a sim- 
ilar model in which the porous matrix was modeled by 
overlapping squares and the tortuosity was computed di- 
rectly as an average over streamline lengths. We found 
that in that model the relation between tortuosity and 
porosity could be approximated by Comiti's and Re- 
naud's logarithmic formula [20] 



T = l-pln(p 



(7) 



with a fitting parameter p ~ 0.77. The dependency of 
the tortuosity on the porosity in the present model of 
overlapping circles is shown in Fig. 7. We fitted the data 
to Eq. (7) and found a good agreement for p w 0.67. 



3.4. Three-dimensional overlapping spheres 

As the final test of Eqs. (5) and (6) we used them 
to find the flow tortuosity in a three-dimensional (3D) 
model of overlapping spheres. The geometry was con- 
structed similarly to the two-dimensional case described 
above. A regular grid of 90 x 90 x 90 nodes was gener- 
ated and freely overlapping spheres of radius 10 l.u. were 
deposited in it to reach the desired porosity value. The 
system was assumed periodic in the x direction and no- 
slip walls were imposed at the four remaining walls. We 
ran the LBM simulation for 10 000 time steps to reach 
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FIGURE 7. Tortuosity as a function of porosity in the two- 
dimensional model of overlapping circles computed using Eq. 
(5). Symbols are our numerical results, the solid line is the best 
fit to Eq. (7). 
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FIGURE 8. Tortuosity, Eq. (5), as a function of porosity 
in a three-dimensional model of overlapping spheres. Open 
symbols are our numerical results, the solid line is the best fit 
to Eq. (8), and the dashed line is the best fit to T = (p^ with 
h fa 0.25. 



the steady state and then computed the tortuosity from 
Eq. (6). 

The tortuosity-porosity dependence in our 3D model 
is shown in Fig. 8 (circles). The solid line in this figure 
depicts a formula that was recently derived analytically 
for a very similar model [21, 22], 
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where B is a constant that depends on the shape of 
the obstacles and the lattice used. The main difference 
between our model and that studied in [21, 22] is that 
we distribute the spheres at random positions and allow 



them to overlap, whereas Eq. (8) was derived for regular 
arrangements of impermeable, non-overlapping objects 
of essentially arbitrary shape. Assuming that Eq. (8) can 
be also used for non-regular arrangements of spheres, we 
fitted our data to this formula using B as a free parameter 
This gave B « 1 .09, which is a bit smaller than B w 1 .209 
derived in [21] for the cubic packing of spheres. As seen 
in Fig. 8, our results are in good agreement with Eq. (8) 
in the whole range of porosities. 
Archie's law [23, 4] 
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(9) 



is another formula for the tortuosity-porosity relation, of- 
ten used for retardation tortuosities, especially the elec- 
trical and diffusional one. Together with Eqs. (2), (3), 
Archie's law suggests T = (p^'^ with b = [n— l)/2. Our 
data can be fitted to this formula with b « 0.25, see Fig. 8, 
which yields « w 1.5, which lies in the range 1 .3 < « < 3 
reported in [4]. A more detailed analysis, involving a 
bigger number of larger systems, is necessary to verify 
which of the equations should be applied to this model, 
and whether Eq. (8) actually applies to non-regular ar- 
rangements of obstacles. 



be expressed as the ratio of the average fluid velocity 
magnitude to the average fluid velocity along the macro- 
scopic flow direction, it turns out to be closely related 
to one of the most fundamental physical phenomena — 
momentum transfer Third, this formula can be applied 
to other forms of transport in porous media, e.g. to dif- 
fusion or electric current [7]. Fourth, our formula can be 
used to flows in the fractal-like [24] or ramified struc- 
tures [25, 26], sytems which are not considered porous. 
For example, one could compute hemodynamical tortu- 
osity in the flow through human artery system, in which a 
fundamental difference between the geometrical and hy- 
draulic tortuosities may be of profound importance for 
medical diagnosis of arterial diseases [27]. 
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4. DISCUSSION AND OUTLOOK 

We have presented several applications of our recently 
introduced method of calculating the tortuosity defined 
as a measure of the average elongation of fluid stream- 
lines [7]. Starting from a simple model of an inclined 
channel, through a U-shaped channel, and ending at com- 
plex 2D and 3D geometries, we found a very good agree- 
ment of this method with other methods serving the same 
purpose. However, the main advantage of our method is 
that it does not require to find any streamlines, which is 
a complicated, time-consuming and error-prone task, es- 
pecially in realistic 3D geometries. Instead, it allows to 
calculate the streamline tortuosity directly from the ve- 
locity field. Not only does this simplify numerical studies 
of this quantity, but should also greatly simplify experi- 
mental measurements of T . 

To summarize, there are several advantages of calcu- 
lating the streamline tortuosity using the method of Ref. 
[7]. First, the formula is simple and flexible, allowing to 
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cal fluid flow system in which the velocity field can be 
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very existence of tortuosity as an average elongation of 
fluid paths. As a consequence, one can concentrate on 
the physical significance of this quantity, including find- 
ing its relation with numerous "tortuosities" defined as 
transport retardation factors. In this context it is inter- 
esting to notice that since the streamline tortuosity can 
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